Here you’ll find my adventures coding pop culture, current events, and maybe if I summon the energy, some helpful lessons from my dissertation work. Right now, I’m a PhD candidate in entomology at Michigan State, but I care about modeling and visualizing topics far beyond insects and plants.
Even though I grew up when Tumblr was in vogue, blogging always scared me until I saw this post and listened to this podcast. Plus, I love sharing my silly little models and visualizations with friends who’ll listen, so why not share them with the world…
Most posts will discuss coding in R, but we may venture elsewhere in the future. Today’s blog post was really fun to put together and was a constant reliving of an amazing day every year—Spotify Wrapped release day!
learn how to use tidytext and apply structural topic models (STM) to identify common trends between songs and years through natural language processing (NLP).
dynamically visualize text and model results with gganimate, plotly, and ggiraph.
spotifyr and the Spotify API
I’ve been a Spotify listener since 2014 and that loyalty paid off for this post with access to seven (!!) years of Spotify Wrapped. If you’re unfamiliar, Spotify releases each user’s top artists and songs for the year every December. Taking the playlists for 2014-2020, I used the spotifyr R package to access the Spotify API and pull out these playlists from library.
Here are some posts that can explain how to do collect data from the Spotify API for non-commerical uses more eloquently than me. The API provides far more information than song names and artists alone; you can obtain audio features and popularity data, and some of the write-ups that dive deep into this topic inspired my post today—tayloR, a spotifyr beginner’s tutorial, & a spotifyr quick rundown.
geniusr and the Genius API
geniusr R package to search the lyrics for each song. This step was tough because some of my top songs are instrumental or remixes and didn’t fit well into my search loop. For those songs, I searched for lyrics separately.
To assemble a dataset that could be useful to interpret with STMs, I had to remove songs were the majority of the lyrics were not in English and filter out stop words like ‘oh’ and ‘ayy.’ A future blog post may model songs with lyrics mostly in Spanish. Since this blog is not about collecting data from the Spotify and Genius APIS, I’m not going to highlight it here. However, to see a skeleton of the code I wrote (minus my playlist IDs and manual song entry), please see this post’s source code here.
tidytext R package to bring my data into forms useful for text analysis. Julia Silge and David Robinson have written a whole book on this, so I’ll direct you to the tidytext book for more information.
I want to take my song lyrics in English and split each song into individual words with unnest_token(). Note: ‘section_name’ refers to chorus, bridge, verse, etc.; ‘row’ is the unique row identifier because some songs appear in multiple years.
library(tidytext)
library(stm)
library(dplyr)
tidy_lyrics_en <- lyrics_df_en %>%
group_by(section_name, section_artist, song_name, song_id, artist_name, row) %>%
mutate(line_number = row_number()) %>%
mutate(song_artist = paste(song_name, artist_name, row, sep = "_")) %>%
unnest_tokens(word, line) %>%
ungroup()
Let’s get rid of most of the words that appear frequently and don’t have much meaning outside the music.
more_stopwords <- c("oh", "ah", "uh", "ooh", "aah", "mu", "woo", "oo", "eh", "di", "da", "ayy", "doo", "doot", "hmm") # I also added some other words to maintain a 'clean' lyric analysis that aren't included here.
lyrics_en <- tidy_lyrics_en %>%
anti_join(get_stopwords(language = "en")) %>%
anti_join(get_stopwords(language = "es")) %>%
filter(word %!in% explicit)
lyrics_en <- as_tibble(lyrics_en) %>%
filter(word %!in% more_stopwords)
To run the structural topic model in the next section, we’ll need our lyrics into a DocumentTermMatrix.
lyrics_dfm_en <- lyrics_en %>%
select(word, song_artist) %>%
count(song_artist, word) %>%
ungroup() %>%
cast_dfm(song_artist, word, n)
To model my Top 100 Songs, we’ll need to determine what is the “best” number of topics (K) that describe the data. These five values for K here are fairly arbitrary (so much as 15 topics is not that different from 16 topics). I spent a long time running a bunch of models not shown here with different values for K.
stms <- data_frame(K = c(5, 15, 25, 35, 45)) %>%
mutate(topic_model = future_map(K, ~stm(lyrics_dfm_en, K = .,
verbose = FALSE, seed = TRUE))) # this code is strongly based on Julia Silge's post
There are several metrics to compare the models varying K including held-out likelihood and semantic coherence. To keep this post from being too long, I only show the residuals to compare the models - but I encourage in something more rigorous to look at all these other metrics.
library(purrr)
library(ggplot2)
k_resids <- stms %>%
mutate(residual = map(topic_model, checkResiduals, lyrics_dfm_en)) %>%
transmute(K,
residual = map_dbl(residual, "dispersion"))
ggplot(k_resids, aes(x = K, y = residual)) +
geom_line()
- It looks like K = 15 topics marginally has the lowest residuals, so we’ll go with that model.
Next, I’ll pull out the values for beta (probability a word is found within a topic) and gamma (probability a topic is found within a ‘document;’ here, the document is a song).
topic_model <- stms %>%
filter(K == 15) %>%
pull(topic_model) %>%
.[[1]]
beta_tidy <- tidy(topic_model)
gamma_tidy <- tidy(topic_model, matrix = "gamma",
document_names = rownames(lyrics_dfm_en))
Here’s the real fun part — interpreting the topics. I love pop music, and pop music doesn’t vary terribly in it’s vocabulary. Nevertheless, I tried my best to label these topics for this heurestic exercise.
CAUTION: I took some liberty “labelling” the topics resulting from the STMs into plain language. There are some quantitative ways to do label topics, like with the labelTopics() function, which I also used to inform my labels. For instance, ‘rain’ and ‘loyalty’ were two of the top words for Topic 6 (labeled ‘longing’) in terms of highest probability and FREX. In Topic 8 (labeled ‘big dreams’), the top words for FREX and Score were ‘stamina,’ ‘dna,’ and ‘greatest.’
library(tidyr)
## Warning: package 'tidyr' was built under R version 3.6.2
topic <- paste0("Topic ", unique(gamma_tidy$topic))
label <- c("party", "confident", "wanderer", "winning", "feeling free", "longing", "airy", "big dreams", "self-love", "friends", "mischievious", "nostalgic", "breaking out", "glow up", "in my feelings")
topic_label <- as.data.frame(cbind(topic, label))
songs <- c("Blue World", "2 Da Moon", "LOYALTY.", "7 rings", "About Today", "Don’t Call Me Up") # we'll dive into these songs later
gamma_abr <- gamma_tidy %>%
mutate(topic = as.character(paste0("Topic ", topic)),
topic = reorder(topic, gamma)) %>%
separate(document, c("song", "artist", "row"), sep = "_", remove = FALSE) %>%
filter(song %in% songs) %>%
distinct(song, topic, gamma) %>%
inner_join(topic_label, by = "topic")
gamma_year <- gamma_tidy %>%
mutate(topic = paste0("Topic ", topic),
topic = reorder(topic, gamma)) %>%
separate(document, c("song", "artist", "row"), sep = "_", remove = FALSE) %>%
inner_join(tracks_full, by = "row") %>%
mutate(artist = artist.x) %>%
select(-artist.x, -artist.y) %>%
inner_join(topic_label, by = "topic") %>%
group_by(year, topic, label) %>%
summarize(mean_gamma_yr = mean(gamma)) %>%
ungroup() %>%
mutate(year = as.numeric(as.character(year)))
## `summarise()` regrouping output by 'year', 'topic' (override with `.groups` argument)
hl_topics <- c("glow up", "wanderer", "in my feelings")
gr_topics <- c("party", "confident", "friends", "wanderer", "good times", "feeling free", "longing", "ethereal", "dreaming", "authentically me", "mischievious", "nostalgia", "breaking out", "in my feelings")
'%!in%' <- function(x,y)!('%in%'(x,y))
gamma_year <- gamma_year %>%
mutate(high = case_when(
label %in% c("glow up") ~ "A",
label %in% c("wanderlust") ~ "B",
label %in% c("in my feelings") ~"C",
label %in% gr_topics ~ "D"
))
First I want a bird’s-eye of my yearly top song data. I wrote this bubble chart resembling the cross-section of a tree trunk:
library(ggiraph)
library(ggplot2)
library(packcircles)
top_artist <- tracks_full %>%
group_by(artist, year) %>%
summarise(count = n()) %>%
arrange(-count) %>%
arrange(year) %>%
mutate(year = as.character(year))
packing <- circleProgressiveLayout(top_artist$count)
dat.gg <- circleLayoutVertices(packing)
colors1 <- as.character(c("#99B898", "#A8E6CE", "#FECEAB", "#FF847C", "#E84A5F", "#D7BDE2" , "#BEC0ED"))
year <- as.character(2014:2020)
colors_year <- as.data.frame(cbind(colors1, year))
top_artist <- top_artist %>%
inner_join(colors_year) %>%
mutate(colors = as.character(colors1))
top_artist$text <- paste(top_artist$artist, "\n", "Song count:", top_artist$count, "\n", "Year:", top_artist$year)
top_artist$id <- c(1:length(top_artist$artist))
dat.gg <- dat.gg %>%
inner_join(top_artist, by = "id")
bubbles <- ggplot() +
geom_polygon_interactive(data = dat.gg, aes(x, y, group = id, fill = factor(id), tooltip = top_artist$text[id], data_id = year), alpha = 0.6, show.legend = FALSE) +
scale_fill_manual(values = top_artist$colors, labels = c('2014', '2015', '2016', '2017', '2018', '2019', '2020')) +
theme_void() +
guides(fill = guide_legend_interactive()) +
labs(title = "Annual rings of my artists' top songs") +
theme(plot.margin=unit(c(0,0,0,0), "cm"),
plot.title = element_text(size = 18, hjust = 0.5)) +
coord_equal()
bubbles_interactive <- ggiraph(ggobj = bubbles, width_svg = 6.2, height_svg = 4.5)
bubbles_interactive
Figure 1. Each bubble and its size represent an artist’s song count in that year. The ‘rings’ of bubbles are colored by year.
A common analysis among the spotifyr posts is to visualize the “mood” of a group of songs with Spotify’s audio features. The “mood” can be understood by the valence (positivity rating where lower values are more “negative”) and energy (lower values have lower energy). If we are to break up a plot of energy x valence like here in the sentify Shiny app, we can see how my songs bunch up into each quadrant.
library(plotly)
library(htmltools)
mood <- tracks_full %>%
plot_ly(x = ~valence, y = ~energy, size = ~rank_weight, sizes = c(1, 350), text = ~artist, hoverinfo = "text")
mood2 <- mood %>%
add_markers(color = ~year, showlegend = F, alpha = 0.05, alpha_stroke = 0.05, colors = colors1) %>%
add_markers(color = ~year, frame = ~year, ids = ~artist, colors = colors1) %>%
add_segments(x = 0.5, xend = 0.5, y = 0, yend = 1, showlegend = F, line = list(color = "#A9A9A9", width = 0.8)) %>%
add_segments(x = 0, xend = 1, y = 0.5, yend = 0.5, showlegend = F, line = list(color = "#A9A9A9", width = 0.8)) %>%
layout(xaxis = list(title = "Valence", showgrid = F, font = list(size = 18)),
yaxis = list(title = "Energy", font = list(size = 18)),
annotations = list(
x = c(.05, .05, .95, .95),
y = c(.05, .95, .05, .95),
text = c("Sad", "Angry", "Pleasant", "Happy"),
showarrow = FALSE,
xanchor = "center",
font = list(size = 15, color = "#547980")
),
width = 600, height = 500
) %>%
animation_opts(1000, easing = "elastic", redraw = FALSE) %>%
animation_button(
x = 1, xanchor = "right", y = 0, yanchor = "bottom"
) %>%
animation_slider(
currentvalue = list(prefix = "Year ", font = list(color = "#696969"))
)
htmltools::div(mood2, align = "center")
Description for topic probability by highlighting some songs.
group_colors1 <- c(`glow up` = "#B6B4C2",
wanderer = "#547980",
`in my feelings` = "#F9CDAD",
`breaking out` = "#AAC9CE",
confident = "#A8E6CE",
`big dreams` = "#CD7672",
`airy` = "#DCEDC2",
`feeling free` = "#63BCC9",
friends = "#E7E7E7",
`winning` = "#7E9680",
longing = "#ECD59F",
mischievious = "#8AB2E5",
nostalgic = "#FFDD94",
party = "#FF8C94",
`self-love` = "#E5C1CD")
ggplot(gamma_abr, aes(x = song, y = gamma, fill = label)) +
geom_bar(stat = "identity", position = "stack") +
scale_fill_manual(values = group_colors1) +
labs(y = "Topic probabilities within song") +
coord_flip() +
theme_classic() +
theme(legend.title = element_blank(),
axis.title.y = element_blank(),
text = element_text(size=12))
Description for topic probability throughout the years
library(gganimate)
group_colors2 <- c(A = "#6C5B7B", B = "#547980", C = "#F9CDAD", D = "#D3D3D3")
group_alpha <- c(A = 1, B = 1, C = 1, D = 0.48)
group_size <- c(A = 1.55, B = 1.55, C = 1.55, D = 1)
time <- ggplot(gamma_year, aes(year, mean_gamma_yr, group = label, color = high, alpha = high, size = high)) +
geom_line() +
geom_segment(aes(xend = 2020.8, yend = mean_gamma_yr), linetype = 2) +
geom_point(size = 2) +
geom_text(aes(x = 2021, label = label), hjust = 0, size = 5.2) +
scale_color_manual(values = group_colors2) +
scale_alpha_manual(values = group_alpha) +
scale_size_manual(values = group_size) +
xlim(2014, 2021) +
transition_reveal(year) +
coord_cartesian(clip = 'off') +
labs(x = 'Year',
y = 'Mean gamma (topic probability)') +
theme_classic() +
theme(legend.position = 'none',
plot.margin = margin(5, 75, 1, 5),
text = element_text(size=20))
animate(time, duration = 9.5)
Reproducible code can be found here (minus the datasets of my Spotify Top Songs and initial cleaning code).